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ABSTRACT 

We present a study of weak lensing shear measurements for simulated galaxy images 
at radio wavelengths. We construct a simulation pipeline into which we can input 
galaxy images of known ellipticity, and with which we then simulate observations with 
cMERLIN and the international LOFAR array. The simulations include the effects of 
the CLEAN algorithm, uv sampling, observing angle, and visibility noise, and produce 
realistic restored images of the galaxies. We apply a shapelet-based shear measurement 
method to these images and test our ability to recover the true source ellipticities. We 
model and deconvolve the effective PSF, and find suitable parameters for CLEAN and 
shapelet decomposition of galaxies. We demonstrate that ellipticities can be measured 
faithfully in these radio simulations, with no evidence of an additive bias and a modest 
(10%) multiplicative bias on the ellipticity measurements. Our simulation pipeline can 
be used to test shear measurement procedures and systematics for the next generation 
of radio telescopes. 
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1 INTRODUCTION 

The bending of light due to the large scale structure in 
the Universe is a powerful probe in studying the underly- 
ing matter distribution. This is due in part to the gravita- 
tional lensing being insensitive to the type of matter causing 
the light deflection (e.g. Bartelmann fc Scfmeider|2001 Re 



frcgicr 2003b, Munshi et al.|2008| >. There is the added benefit 
of lensing being sensitive to the geometry of the Universe, 
making it useful in the study of dark energy ( Hutererf2 002). 

To date, almost all cosmic shear analyses have been 
conducted at optical wavelengths. However, radio astron- 
omy is currently going through a period of rapid expansion 
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which will make future radio surveys competitive for lensing 
studies. New radio telescopes will have sensitivities that will 
reach a level where the radio emission from ordinary galaxies 
will be routinely resolved (e.g. with eMERLIIs Q LOFAB0 
and eventually SKA[^] c.f. Seymour et al.|2004| ; hence radio 
source densities will become comparable to those found at 
optical wavelengths. Also, radio interferometers have a well 
known and deterministic dirty beam pattern, which may be 
an advantage in deconvolving galaxy shapes for shear mea- 
surements. 

With the FIRST survey ( |Becker et al.|1995|[White et al 



1 http://www.e-merlin.ac.uk/ 

2 http://www.lofar.org/ 

3 http://www.skatelescope.org/ 
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19971, Chang et al. (2004) made the first detection of cosmic 



shear with radio data. This survey has a detection thresh- 
old of 1 mJy, with ~ 20 resolved sources per square degree 
useable for weak lensing. This is a much lower number den- 
sity than found in deep optical shear surveys, with a cor- 
respondingly lower signal-to-noise on the final cosmological 
constraints. However, the differential radio source counts at 
1.4 GHz show an increase at flux densities below 1 mJy, (e.g. 
Seymour et al.|[2004 ), and it is this increase in the number 
density at the micro-Jansky level that makes future radio 
weak lensing plausible. 

The feasibility of weak lensing studies at radio wave- 
lengths, and in particular at the micro-Jansky flux levels, 
was demonstrated in Patel et al. (20101. Due to the low 



number density of sources used in that work, no significant 
cosmic signal was detected. However, an upper bound was 
obtained on a combination of the cosmological parameters 
erg, the normalisation of the matter power spectrum, Q m , 
the cosmic matter density parameter, and z m , the median 
redshift of the sources. One of the main conclusions of that 
work was the need for a detailed study of the systematics 
involved in radio interferometry and the relevant imaging 
techniques. 

In this current paper, we pursue that study. We will ex- 
plore possible systematics which will be important for weak 
lensing studies with future radio surveys; for example, the 
systematics which might be introduced by steps such as the 
CLEAN algorithm or a poor choice of shape measurement 
parameters. In tandem, we will quantify our current abil- 
ity to measure realistic galaxy shapes from simulated radio 
data, increasing our confidence in current and future radio 
lensing measurements. 

The paper is organised as follows: in ^2]we describe the 
whole of our radio image simulation and shape measure- 
ment pipeline. We describe the shapelets method which is 
used throughout this work, and how deconvolution works 
within its framework and the shape estimator we use. We 
then describe the shapelets based images that we created as 
the input for our simulations, before describing how the sim- 
ulator works to produce realistic radio interferometer images 
that we use for our analysis. We also describe the telescope 
configurations we use in this study, and the observational 
effects that we consider. 

In S|3]we describe the results of our shape measurements 
on the simulated images. We assess the appropriate level of 
CLEANing of images, and examine the modelling of point 
sources and the effect of changing the position on the sky 
at which images are observed. We discuss the impact of 
the shapelet scale parameter in modelling galaxy elliptici- 
ties successfully, and the effect of slightly changing the uv 
sampling. We then present the main result of the paper: 
how well we are able to recover the input ellipticities with 
our shape-measurement method. We describe how we add 
realistic noise into the simulations and what effect this has 
on our results. In ^4]we end with a discussion of our results, 
their limitations, and suitable directions for further study. 



2 SIMULATIONS OF RADIO IMAGES 

The radio image simulation pipeline is built in two parts. 
We first create a suite of input images, containing sources 



which constitute the inputs for the main part of the simula- 
tion code. These input images are then fed into the simulator 
that 'observes' and 'images' them with a given telescope con- 
figuration and user defined observation parameters. In this 
section we start by briefly describing the shapelets technique 



of Refregier ( 2003a I and Refregier & Bacon ( 2003 1 which we 



use to create the galaxy shapes in our input images; we will 
also use shapelet techniques for shape measurement later. 
We then describe the creation of the input images before 
discussing the radio observation pipeline. 

2.1 Shapelets 

The shapelets description of galaxy shapes is based on de- 
composing an object's surface brightness /(x) into a series 
of localised basis functions -B nil „ 2 , called shapelets: 



/( X ) = /»1,"2 B «1,» 2 ( X ^)- 



(1) 



We briefly describe the relevant parts of the method here 
and refer the reader to Refregier I ( |2003a[ ), [Refregier fc Bacon 
H2003D and |Massey fc Refregier] ( |2005[ ) for a more detailed 
description. The basis functions used (in the Cartesian for- 
malism) B nit „ 2 , are a localised orthonormal basis set: 



B„ 



.(x;0) = 



H ni (f)H n2 (f) 



(2) 



(2 n i n 2/3 2 7rni!n 2 !) 2 

where H m (f3) is the m th order Hermite polynomial, and /3 
describes the characteristic scale of the basis. We will see 
that this quantity is very important in our PSF and galaxy 
shape analysis, ni + ri2 refers to the order of the basis func- 
tions, and in practice all galaxy decompositions are trun- 
cated at some order n ma x = n\ + n^. n max needs to be 
chosen so that the galaxy model is sufficiently detailed to 
capture ellipticity information. As the basis functions are 
orthogonal, we can find shapelet coefficients for a galaxy by 
calculating 



/(x)B ni ,„ 2 (x;/3)d x. 



(3) 



Massey & Refregier ( 2005 ) introduced the closely related po- 



lar shapelet formalism (the polar basis set is an orthogonal 
transformation of the Cartesian set; see Massey & Refregier 
(2005) for details of this transformation), which has a dif- 
ferent basis set P n ,m,(0,4>> P), given by 



P n , m {9,<t>;P) 



(-1) 



(n-|m|)/2 



3 • 



-M 



[(n-H)/2]! 
T[(n-|m|)/2]l 



e -0'/2/3* e -im<t> 



1/2 



(4) 



(n-\m\)/2 ' 

In this basis set 6 is the modulus of the complex sky posi- 
tion vector 9\ + i02, and <j> = arctan(#2/#i). Li(x) are the 
associated Laguerre polynomial defined as 

x-"e x <¥ 



L*(x) 



Two integers, n and m uniquely describe every member of 
the basis set, with n > and \m\ ^ n. The surface brightness 
of a galaxy f(9) is given by 

oo n 

f(e)=f(9,4>)=J2 E fn,mPn,m(0,4>;P) (6) 
n— m— — n 



') 



(5) 
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The polar basis set has rotational symmetries which are very 
useful for describing weak lensing, so in practice galaxies 
are decomposed using the Cartesian basis set (which easily 
describes square pixels) and then transformed to the polar 
set. 

Both the Cartesian and polar shapelet basis functions 
have simple behaviour under convolution (Rcfrcgicr 2003a I 
and deconvolution ( jRefregier fc Ba con 2003), making them 
particularly well suited for describing and correcting the ef- 
fects of a PSF. Since this is important for our shape mea- 
surement analysis, we briefly describe deconvolution in the 
shapelet framework below. 



2.2 Deconvolution with Shapelets 

The approach used in this work for deconvolution is to es- 
timate the deconvolved shapelet coefficients f ni ,n 2 by 'for- 
ward convolving' the shapelet basis function with the PSF 
model g(x), in advance, creating a new basis set which we 
label 



t(x;P) = 3(x;/3) * B„ 1 ,„ 2 (x; / 8), 



(7) 



with an equivalent expression for the polar shapelet basis 
functions. Fitting the data ft(x) using this new basis set 
D ni ,n 2 , on e obtains the deconvolved shapelet model for the 
galaxy as follows: 

/i(x) = £?(x) * /(x) 

OO 

= y( x )* Yl f^un 2 B nijn2 (x;fi) 

OO 

= /fti,n a bW*5 ni)fta (x; i 8)] 

Til , n 2 
OO 

= ^ /n 1 ,n 2 -Dn 1 ,n 2 (x;/3) (8) 

Comparing with Equation [T] the returned coefficients will 
reconstruct the the deconvolved image when used with the 
original basis set. 

We use the p ublicly available shapelets so ftware pack- 
ag^] described in Massey & Refregier ( 2005 1 in order to 
make shapelet decompositions for all our objects. This code 



is well tested using optical data (c.f. Heymans et al. 2006 



and Massey et al.|2007a[). In Patel et al. ( 2010 \ we indicated 



the applicability of this code to radio data with some modi- 
fications; with our radio image simulations, we are now in a 
position to demonstrate its ability to accurately recover el- 
lipticities in i|3] In all that follows we have used the shapelet 
equivalent to Gaussian weighted quadrupole moments to es- 
timate ellipticities for all our objects. In the polar shapelet 
convention this estimator takes the form 



£1,2 



{fo,0 - /4,o) ' 

as is fully discussed in |Massey et al.| |2007b[ ). 

4 http:/ /www. astro. caltech.edu/~rjm/shapelets/code/ 



(9) 



2.3 Image Simulations with Shapelets 

In this section we describe how the input images are created. 
We generate two different sets of input images: the first set 
contains point sources, consisting of a collection of single 
illuminated pixels at random locations. The second suite are 
images containing detailed galaxy shapes and morphologies, 
but each with a centroid at the same locations as the point 
sources in the previous image set. We will use the point 
sources to study the behaviour of the dirty beam or point 
spread function of the telescope. The galaxy images provide 
a more realistic challenge, with the goal being to recover 
the (known) shapes of galaxies in the presence of all the 
systematics that distort them. 

The galaxy images are created using the shapelets based 
formalism described above. The method of populating the 



images is based on the Massey et al. ( 2004J) pipeline, which 
was used as part of the STEP2 shear methods testing pro- 
gramme ( |Massey~et al. 2007a). The motivation for the use 
of shapelets galaxy models, as discussed in |Massey et al.| 
(2004), is to simulate deep sky images which include some 
degree of realistic galaxy morphology. 

The procedure used to generate simulated galaxy im- 
ages in this paper most closely follows that described by 



|Rowe et al. ( 2012 1 for the generation of simulations of Hub- 
ble Space Telescope (HST) data. Using a PSF and real noise 
properties estimated directly from real HST survey data 



(specifically the GEMS Survey: see, e.g., Heymans et al. 
[2505} |Rix et al.|[2004| , |Rowe et al^ $20lty created realistic 
simulations of Advanced Camera for Surveys (ACS) images 
in the F606W filter. These were then used to investigate 
shear and higher order lensing shape measurements in ACS 
data. 

As we are simulating radio observations, we do not re- 
quire a high level of similarity to optical images, in terms 
of noise, PSF or even shape; in our previous work, Patel 
et al. (20101, we found that on a galaxy- by- galaxy basis, 



optical and radio shapes of galaxies are only weakly corre- 
lated. However, we also showed that the overall distribution 
of ellipticities in the optical and radio sky were very similar, 
so we choose to follow the optical distribution of shapelets 
in this study. This is a reasonable choice given the current 
absence of a large observed sample of highly resolved radio 
objects at fijy flux density thresholds, which could inform 
us about the details of the radio shapelet distribution. 

A 'starter set' of galaxy shapelet models was created 
using shapelet decompositions of GEMS images, modelled 



according to Rowe et al. (20121 using a spatially and tem- 



porally varying model of the ACS PSF. The temporal vari- 
ation of each shapelet coefficient in an n max = 20 model of 
the PSF was modelled in three epochs as described by |Hey-| 
|mans et al. (20051, using a third order polynomial surface 
to describe the spatial variation on each ACS chip. The re- 
sulting PSF-deconvolved shapelet models are then sampled 
from to produce the simulated galaxy images. 

When creating the simulated galaxies from the starter 
set, the models are randomly rotated and/or inverted to 
eradicate any remaining signature of gravitational lensing. 
We also note that the starter set represents a large but ul- 
timately limited sample of galaxy morphologies. This is al- 
leviated by introducing small random perturbations to the 
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shapelet models (see Massey et al.|2004 & Rowe et al.|2012 
for details). 

Using this prescription, a suite of 100 images of size 
1024 x 1024 pixels with a pixel scale of 0.05" were created. 
We have constructed these images such that all the sources 
(galaxies and point sources) were at a constant flux. Galaxies 
have been assigned the same flux, as the process of calibra- 
tion should in principle remove all brighter galaxies in a ra- 
dio image. The brighter galaxies furthermore will be present 
in smaller numbers per square arcmin, and hence produce 
less side lobe noise. Therefore it is the galaxies which are 
at the sensitivity limit of the survey which will dominate 
the side lobe contribution to the image. We could have pro- 
duced a simulation with a luminosity function, but this is 
unknown at the faintest magnitude limits, so the approxi- 
mation we have made should suffice. 

We also inspected the relationship between the /3 values 
of the simulated galaxies, derived from GEMS, and the f3 pa- 
rameters that were derived from the shapelet modelling of 



the actual radio sources in Patel et al. (20101. It was found 



that the distribution of the /3 parameters from the best- 
fitting shapelet models of the GEMS data had a systemat- 
ically lower value (shift of 0.185") than the radio shapelet 

0.3" for the silver 



models from 



Patel et al. 



(2010| ((/?) 

catalogue in that work). To make the simulations a more re- 
alistic representation of radio data we therefore add 0.185" 
to the /3 scale parameter for our galaxy models; these then 
provide a close match to the distribution of the derived /3 



values from Patel et al. (2010) 



For each image, we have a corresponding catalogue con- 
taining the input shapelet coefficients of each source in the 
image, as well as the centroid of the objects. An example of 
one of our images containing galaxies is shown in Figure [I] 
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senting deep pointings. The input images contain a number 
density of n ~ 140 arcmin -2 , which is far larger than current 
lensing studies, even those in space. The reason for setting 
this high number density is to probe systematic errors re- 
lated to radio interferometers deriving from both the instru- 
mental setup and the radio imaging pipeline. While a large 
number of point sources allows one to better probe the dirty 
beam behaviour across the field, a larger number density 
in the galaxies means a more challenging imaging problem, 
since sidelobe noise will be greater in densely packed images. 

Note that the simulated galaxies have not been sheared; 
in this analysis we are not concerned with trying to recover a 
cosmic shear signal from galaxies; rather, we are examining 
how well we can recover known input ellipticities, which is 
the basis for recovering such a signal. 

Having described the details of the input images, we 
now describe the details of our simulated configurations. We 
aim to create observations from radio arrays of our choice, 
including imaging steps such as dirty beam deconvolution 
and visibility weighting, in order to create realistic restored 
images. 



2.4 Simulations Pipeline 

Our simulation pipeline is implemented using the MeqTree^] 
software system. MeqTrees is a software package for imple- 
menting the Measurement Equation for radio interferome- 



ters and is fully described in Noordam & Smirnov (20101 
and |Smirnov| ( |2011| > . The heart of a Measurement Equa- 
tion is formed by the 2x2 Jones matrices (Jones 19411 



which describes the various effects associated with obser- 
vations that can corrupt the measured visibilities. The for- 
mulation for a generic Radio Interferometer Measurement 



Equation (RIME) was developed by Hamaker et al. (19961, 
after preparatory work by |Morris et al.| ( |1964[ ). |Hamaker| 
(20001 then recast the formalism into 2x2 matrix form 
which is used within MeqTree^] The RIME provides an ele- 
gant mathematical framework for generic radio instruments, 
both existing and future, to be better understood, simulated 
and calibrated. For a full description of the mathematical 
formalism for the measurement equation of interferometers 



we refer the reader to Hamaker et al. ( 1996 1, Smirnov ( 2011 1 



Figure 1. Example of input shapelet-based galaxy images; the 
image has a linear greyscale. 



and references therein. 

The MeqTrees software was originally designed to im- 
plement the measurement equations for the purposes of sim- 
ulation and calibrating ( Noordam fc Smirnov|2010 1 . We use 
MeqTrees in this work to simulate the sky as would be ob- 
served by a specific radio interferometer. The galaxy images 
created above are fed into the MeqTrees simulator, which 
calculates observed visibilities; these visibilities are then im- 
aged. The user can specify many options such as the level of 
noise on the visibilities, the specific method used to decon- 
volve the dirty beam, and the weighting scheme applied to 
the visibility data prior to imaging. 

One subtlety is that the input in our work is an image 
rather than a set of visibilities. In order to obtain uv samples 
from the image, we use a module of the MeqTrees software 



We note that these images represent an ambitious imag- 
ing scenario, as they are very densely populated fields repre- 



5 http://www.astron.nl/meqwiki 

6 Some versions of the H1ME are still implemented using the 4x4 



■ 



Mueller matrices ( Mueller [1948 I , which are entirely equivalent 
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Figure 2. Summary of the main steps involved in the simulation pipeline from start to finish. 



that implements a 'de-gridding' algorithm. This algorithm 
is an interpolation scheme that allows one to transform be- 
tween the regularly gridded image and a sparsely sampled 
Fourier (or uv) plane. The details of how this algorithm work 
is explained fully in (Tan 1986). 



The well known relationship between observed visibili- 
ties and the true sky brightness for an interferometer is given 
by: 



I(£, m) 



V(u, v)e 



2i7r(u£-\-vm) 



dudv, 



(10) 



Equation |10| should also include the sampling function term, 
S(u, v) on the left hand side. The sampling function ac- 
counts for the fact that the visibilities are collected at a 
set of discrete locations in the uv plane. Its functional form 
is simply a linear combination of Dirac delta functions at 
all (it, v) where data is collected. The Fourier inversion of 
Equation [To] leads to the dirty image, which is the convolu- 
tion between the sampling function (in the image plane) and 
true sky brightness. Since Direct Fourier Transforms become 
computationally expensive with large datasets, Fast Fourier 
Transforms (FFTs) need to be used. Although the use of 
FFTs speeds up computations, it requires that we grid the 
data. The process by which this de-gridding is done, and 



its implication for the recovered image is described in |Tan| 
(19861. The MeqTrees software employs the use of prolate 



spheroidal functions in the de-gridding process, which have 
been shown to reduce the effect of aliasing ( Tan|1986 i. 

In Figure[2]we illustrate the simulation setup from start 
to finish. (Note that in this work we have focused on one par- 
ticular method of deconvolution of the dirty beam, namely 
the CLEAN algor ithm ( |Hogbom|[l974[ ), which we shall de- 
scribe in Section 2.5.21. We start with an image that is 



Fourier transformed to create the set of simulated visibil- 
ities. These are then modified by the RIME formalism for 
our given set of parameters to create a set of visibilities as 
observed by our given array. These visibilities are Fourier 
transformed again to produce the dirty image, which is then 
CLEANed as our setup defines to produce the CLEAN com- 
ponent image. The residual of the dirty image is then added 
to the CLEAN components once convolved with the CLEAN 
beam (usually an elliptical Gaussian measured from the 
main lobe of the dirty beam). The final step involves con- 
volving the CLEANed image with the CLEAN beam. The 
CLEAN beam is usually a Gaussian fitted to the main lobe 



of the PSF and it is this final beam that the point source 
simulations will characterise. This is also what we shall de- 
convolve from the galaxies to produce our final images from 
which the shapes of our sources will be measured. 



2.5 Simulated Configurations 

In this section we describe the different experiments we have 
carried out to examine a variety of potential causes of sys- 
tematic effects. 

We choose to concentrate on two specific interferome- 
ters: eMERLIN and LOFAR. The eMERLIN array is an up- 
grade to the existing MERLIN array, maintaining the same 
number of dishes. The eMERLIN array upgrade is designed 
to increase sensitivity by more than an order of magnitude 
by using new receivers and telescope electronics. The array 
spans 217 km, with observing bands at 1.3-1.8 GHz, 4-8 GHz 
and 22-24 GHz, with a total bandwidth of 4 GHz. It also has 
resolution capabilities of between 10-150 mas and sensitivity 
of ~ 1/^Jy. 

LOFAR is a new generation radio interferometer, with 
the ultimate goal of surveying the Universe at low frequen- 
cies with high resolution and sensitivity. It operates in the 
less explored low frequency range of 10-250 MHz. LOFAR 
belongs to a new generation of telescopes, with the concept 
of utilising many inexpensive dipole antennae arranged in 
stations without any moving parts, in contrast to the usual 
notion of radio dishes as used by eMERLIN. Different parts 
of the sky are observed by steering the beam electronically. 
The spatial resolution of the Netherlands part of the array is 
governed by the ~ 100 km baselines, leading to a resolution 
of a few arc minutes at 240 MHz. The LOFAR array has 
now been extended over Europe, with stations in the UK, 
Germany, Sweden and France. These Europe wide baselines 
reach ~ 1500km, leading to sub-arcsecond resolution above 
100 MHz. In our study we include these larger, Europe- wide 
baselines. 

For both arrays we create two different Measurement 
Set|](MSs) centred on a fixed RA 00 h 02 m 34 s .43, but with 
differing DEC +90°16'4l".75 and DEC +60°16'4l".75. A 



Interested 



readers 



referred 



to 



|http:/ /aips2.nrao.edu/docs/not es/229/229. html| for a detailed 
explanation. 
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LOFAR uv COVERAGE, <5=+60* 



LOFAR uv COVERAGE, 6= +90° 
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Figure 3. uv coverage corresponding to our 4 starting MSs. The top panels show LOFAR coverage for observations at two declinations, 
while the bottom panels show eMERLIN coverage. 



Measurement Set is a specific definition of how visibility 
data is stored, designed to be compatible with the measure- 
ment equation formalism. For our purposes we can think of 
it as a large table containing information about the particu- 
lar observation in question. For the eMERLIN MSs we use a 
central frequency v = 1.4 GHz while for the corresponding 
LOFAR MSs we use a central frequency v = 240 MHz. In 



both cases we set our bandwidth specifications to 128 chan- 
nels each of 125 kHz resulting in a 16 MHz bandwidth. We 
have also employed a 20 sec time averaging over a 24 hour 
observation. The uv coverage corresponding to these four 
MSs are shown in Figure [3] 

We use the MeqTrees imager to create restored (decon- 
volved) images for these telescope configurations. We image 



Weak Lensing Measurements in Simulations of Radio Images 7 



an area twice the size of the original input image, and use a 
portion of this extra image to estimate the rms of the residu- 
als as described below. Our default imaging option uses the 
Clark CLEAN deconvolution algorithm; we shall describe 
below how we estimate the number of CLEAN components 
used for each simulation. In creating the restored images we 
have kept our pixel scale the same as that of the original 
input images, 0.05". 

In some runs of the analysis we do not include mea- 
surement noise (i.e. the visibilities are not corrupted in any 
way). Clearly this is unrealistic, but allows us to distinguish 
between systematic effects and noise effects. Once we have 
analysed these noise-free images, we explore the effect of 
adding measurement noise to see what effect that has on 
our results. Below we list and briefly discuss the telescope 
effects we study in this work. 



2.5.1 Observing Angle 

We have created MSs for both arrays observing at differing 
declination angles to examine what effect this has on our 
shape measurements. It can be shown that the equation for 
the ellipse traced in the uv plane by a particular baseline is 
given by ( |Thompson|T999l ) 



that 



I D (£,m) = ^2AiB(l- 



mi)+h{l,m), (12) 



2 , 
U + 



v — (Lz /A) cos So 
sin So 



L'x + L'y 
A 2 



(11) 



where Lx,y,z are the coordinate separations between the 
two antennae and So is the declination of the phase tracking 
centre (usually where the field-of-view is centred). As the 
interferometer observes a point on the celestial sphere, the 
rotation of the Earth causes the u and v components of the 
baseline to trace out an elliptical locus according to this 
equation. 

Since the sampling function is effectively a collection of 
these elliptical loci (c.f. Figure I5J, we see that it depends on 
the declination of the observation and the antenna spacings. 
The sampling function, as already mentioned, determines 
the PSF (or dirty beam) for the experiment. Our different 
declination observations give us the means to see how the 
variation of the sampling function (measured from the point 
sources) affects the recovered galaxy shape distributions. 



2.5.2 CLEAN 

We have mentioned above the Fourier relationship between 
the observed visibilities and the desired sky intensity. We 
discussed how the FFT of the visibilities leads to an esti- 
mate of the dirty image, which is the sky brightness con- 
volved with the sampling function. In order to obtain an 
estimate of the sky intensity we need some way to perform 
a deconvolution. 

The commonly used algorithm CLEAN is one way to 



perform this deconvolution. Devised by Hogbom ( 1974 1 it 



assumes that the radio sky can be represented by a collec- 
tion of point sources (CLEAN components) in an otherwise 
empty field. The intensity distribution 1(1, m) is then ap- 
proximated by superposition of these point sources which 
have a positive intensity Ai, at the locations (£i,mi). The 
CLEAN algorithm then aims to determine Ai(li,mi) such 



where I D (£, m) is the dirty image that is obtained from the 
inversion of the visibilities, B(£,m) is the dirty beam which 
is the inverted sampling function and, h (£, m) is the resid- 
ual brightness distribution. The approximation is deemed to 
have been successful if this residual noise is similar to that 
of the measured visibilities. This decomposition cannot be 
analytically computed and an iterative approach is required. 

The original CLEAN algorithm is applied entirely in 
the image plane. It uses a simple iterative procedure to find 
the position and strengths of the sources in the dirty image, 
from which a dirty beam multiplied by the peak strength 
is subtracted. All positions where this occurs are recorded 
as well as the corresponding peak flux. The procedure stops 
when all remaining peaks fall below some specified level. 
The recorded positions and fluxes constitute the point source 
model, which is then convolved with the idealised CLEAN 
beam (usually the central lobe of the dirty beam) and the 
residuals from the dirty image are added to produce the 
final, deconvolved image. We refer the reader to |Schwarz| 
|1978| and |Schwarz1 ( [Tff9| for further details. 

Our imaging pipeline can implement a variety of 
CLEAN algorithms; in this paper we choose a widely used 
version, the Clark CLEAN (|Clark||1980|) , which efficiently 



implements the algorithm and can provide improved speed 
for larger images. 

It is important for us to study the effect of the CLEAN 
method (and the process of deconvolution as a whole) on 
shape measurement. Although CLEAN is well tested and 
appropriate for many imaging applications, it is non-linear 
and does not necessarily converge in a well-defined manner. 
It is not immediately clear if it is suitable for the purpose 
of shape measurement as required for weak lensing; we will 
test this in the next section. 



2.5.3 Sampling of the uv Plane 

In creating the MSs, we have examined 2 different sampling 
configurations. As described above, our default MSs were 
created using a regularly sampled uv coverage by sampling 
a 24 hour observation every 20 seconds. Further to this we 
created an irregularly sampled uv coverage by sampling a 
1000 hour observation every 20 minutes. This provides a 
similar density of visibility samples to the previous config- 
uration, but these samples are not uniformly spaced along 
the tracks. Although this particular configuration is unlikely 
to arise in a normal observation, irregular sampling of the 
uv plane can come from, for example, data flagging at some 
given interval. Since the uv coverage determines the PSF, 
which we know to be important for weak lensing studies, it 
is useful to see if the results are sensitive to the sampling 
configuration. 



3 SHAPE MEASUREMENT RESULTS 

In this section we present the results from the shape mea- 
surement analyses of both point sources and galaxies. We 
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firstly address the question of how many CLEAN compo- 
nents are required to adequately represent these images. 
We then examine the behaviour of point source elliptici- 
ties, for different telescope configurations and declinations. 
Next we assess the modelling of galaxy shapes, in relation 
to the shapelet scale parameter and uv sampling. All this 
is in preparation for the main result of the paper, which is 
the presentation of whether output ellipticity estimates are 
in line with true input ellipticities. Finally, we discuss the 
impact of realistic noise on the measurement of ellipticity. 

3.1 Required Number of CLEAN Components 

We aim to find the lowest number of components that suit- 
ably deconvolve the sources in the images, since this reduces 
computation time and avoids fitting noise. To assess this 
issue, we set several different target numbers of CLEAN 
components on one of the images from our set of simula- 
tions. We ran this experiment on all four MSs with the point 
sources and the galaxies. Since the point sources have a much 
simpler set of source structures, it would be expected that 
they would require fewer CLEAN iterations. We examine 
the residuals (of the dirty image) in order to compare the 
relative merits of the different numbers of CLEAN compo- 
nents used. We examine where we find the minimum rms 
residuals as a function of number of CLEAN iterations. 

We image the residuals just outside the patch where we 
have sources in order to estimate the residual side lobe noise 
from the sources inside the central part of the image. To es- 
timate the rms of the residuals, we select a frame around the 
inner image containing the sources, with a width 16 pixels 
(0.8") which we split into 20 cells, from which we calculate a 
mean rms. In Figure [4] we show how the mean point source 
and galaxy rms residuals vary with the number of CLEAN 
components for each of our MSs. 

We wish to perform our analyses on the images with the 
smallest rms residuals. For the point source images, the low- 
est tested number of CLEAN components, Ncln = 10 5 , has 
smallest rms, and we adopt this value for point source analy- 
sis. We see that in the eMERLIN galaxy simulations there is 
a minimum in the rms residuals near Ncln = 3 x 10 s , while 
for the LOFAR cases the minimum lies near Ncln = 2x 10 , 
these are the iteration numbers that we have adopted for 
galaxies in the rest of this paper. We note that these values 
are not dependent on the declination of the observations. 

3.2 Source Extraction and Shapelet Modelling 

Since the images we create have known source positions, 
we bypass the source extraction stage of a lensing analysis. 
We instead use SExtractor in the ASSOC mode whereby we 
input a catalogue of positions where objects are known to 
lie. Of course, even in this mode we will only detect objects 
which meet the internal requirements to be classified as an 
object. We use the default SExtractor settings of Table [l] 
with these settings we find that we are able to detect ~ 
85 — 90% of the objects that are in the original images, 
with detected number densities ranging from n ~ 120 — 130 



Table 1. Tabic summarising main SExtractor parameters. 



Configuration Parameter 



Value 



ASSOC.TYPE 
ASSOCJtADIUS 
DETECT.THRESH 
DETECTJV1INAREA 
DEBLEND_MINCONT 
DEBLEND_NTHRESH 
BACK_SIZE 
BACK_FILTERSIZE 



NEAREST 
5 pixels 
1.5 

5 

0.005 
32 
61 
3 



the beam, which we can then use to deconvolve from the 
galaxies. After this, the deconvolved galaxy ellipticities will 
be compared to the ellipticities of the original catalogues. 



3.3 Point Sources 

For the point sources we use a straightforward approach for 
shapelet modelling, decomposing each source into shapelets 
while examining how the distribution of ellipticities changes 
if we vary the /3 parameter. We try the algorithm employed 



We use the shapelets software as described in §2.1| to 
estimate the ellipticities of the point sources and the galax- 
ies. We first examine the point sources, in order to estimate 



in the shapelet software ( Massey & Refrcgicr 2005 1 to op- 
timise the ft and n max parameters, as well as adopting /3 
values of 1, 1.5 and 2 pixels with a fixed n ma x = 10. We find 
the mean point source FWHM is approximately ~ 2 pixels 
for the eMERLIN simulations and ~ 6 pixels for the LOFAR 
ones (where the pixel size is chosen to be the same as for the 
input images). In Figure|5]we show the normalised ellipticity 
distributions of the point sources for each measurement set; 
note that the 'Massey' label refers to ellipticities obtained 
using the optimisation algorithm in the shapelets software 
I Massey fc Refregier|2005 1. 

From these ellipticity distributions we can see that the 
most tightly peaked distributions are produced from a fixed 
/3 = 1.5 pixel approach. The underlying PSF ellipticity dis- 
tribution should be narrow in a small image region, so we 
adopt P = 1.5 for our PSF measurements in all four cases. 
We also note that in the LOFAR ellipticity distributions, 
changing /3 can change the measured ellipticity of the PSF 
and therefore would lead to different galaxy ellipticities if 
used for deconvolution. This can be understood by com- 
paring the mean point source FWHM of 6 pixels with the 
attempted /3 choices; having a small /3 will mean that the 
decomposition cannot capture the larger scale shape infor- 
mation, so the /3 — 1 pixel histograms behave erratically. 
Since the eMERLIN mean point source FWHM is ~ 2 pix- 
els, all our attempted (3 choices produce well measured point 
source models. 

To quantify any spatial variation in PSF ellipticity, we 
combine all the point source simulations and bin objects' 
ellipticities in x and y coordinates. In Figure [6] we show 
how the mean ellipticities vary as a function of pixel posi- 
tion across our field. For each of our four MSs we show the 
measured (uncorrected) and corrected ellipticities; the cor- 
rection applied in each case is a simple subtraction of the 
global mean measured ellipticity on the image. 

In all four cases, we see that there is no substantial 
variation in the beam behaviour across the images. In each 
case we find that e 2 ^ 0. For the LOFAR DEC = +90° 
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Figure 4. Mean rms residuals for varying number of Clark CLEAN components, for point sources and galaxies, in each of our four MSs. 



there is a small ei component with mean 0.005, while in 
the lower declination LOFAR simulation we see there is an 
almost constant ei = 0.06 across the field. After correction, 
there is no evidence of a coherent ellipticity pattern. For 
the eMERLIN observations, we see a change in sign of the 
magnitude of ei from positive to negative for the DEC = 
+90° and DEC = +60° cases respectively. As in the LOFAR 
case, the induced ellipticities are constant in our small fields, 
so a simple correction by the mean works well. 



Table 2. Table summarising PSF ellipticities. 



MS 




PSF Ellipticities 






ei 




LOFAR <5 = 


+90° 


0.005 


-10~ 4 


LOFAR <5 = 


+60° 


0.063 


-10~ 6 


eMERLIN 8 -- 


= +90° 


0.024 


-0.001 


eMERLIN 5 -- 


= +60° 


-0.042 


-10" 5 



We can now use these models from the point sources 
to construct a PSF in each of our four measurement sets 
and use them to deconvolve the galaxies. The peaks in the 
histograms provide us with a tight measure of the beam el- 
lipticity, as does using the mean of the uncorrected whisker 
plots; the measured ellipticities for the beam are summarised 
m Table [2] for our four MSs. To construct our PSF models, 
we stack all our point sources and perform a shapelet decom- 
position with n max — 10 and /3 = 1.5 pixels on the stacked 
source. 



3.4 Shapelet Modelling the Galaxies 



In Patel et al. (20101 we found that the optimisation algo- 
rithm from |Massey fc Rcfrcgier (2005) performed poorly on 
galaxies in our radio data, and that fixing the n ma x and ft 
parameters in the shapelet modelling gave much better re- 
constructions. We attributed the poor performance of this 
algorithm to the properties of the noise in the radio data. 
Here, we have further tested this by performing the shapelet 
decompositions of galaxies in a number of different ways: we 
have used the optimisation algorithm from Massey & Re- 
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Figure 5. Normalised ellipticity distributions for our four default measurement sets, for different methods of shapelet decompositions. 
Red curves correspond to the shapelet optimisation algorithm of (Masscy & Rcfrcgicr 2005), and the yellow, blue and green curves 
correspond to decompositions done with fixed n ma x = 10 and (3 = 1, 1.5, 2 pixels respectively. The solid lines show ellipticity element ei 
and the dashed show t2- 



fregier ( 2005 1 to see if it performs better on the noiseless 
simulated images; we have also performed several shapelet 
decompositions with a fixed n max = 10 and varying /3; we 
have used multiples of the SExtractor FWHM for this pur- 



pose, analogously to Patel et al. (20101 



As with the point sources, we use SExtractor in ASSOC 
mode for source extraction before performing the shapelet 
decompositions. Before comparing to the original input el- 
lipticities we removed all shapelet model failures. Most of 
these shapelet models can be attributed to objects being 
close to other objects and also objects lying near the edge 
of the images. Due to the high number density of objects in 
the images there is a large failure rate; in all four cases we 
lose ~ 50% of objects due to bad shapelet models. 

For all our catalogues we compute the deconvolved el- 
lipticity estimator described in §2.1[ with a PSF kernel with 
P = 1.5 pixels. We then match our catalogues to the orig- 
inal catalogues and compute the corresponding 'true' ellip- 
ticities, e| to which we will compare. We have binned the 
catalogues in e* and have calculated the median measured 
a in each bin, along with associated error. 



We can now assess the best choice for the /3 parameter 
for galaxy shapelet decompositions. For each choice of /3 we 
fit a linear model to our data points, (e.i — e\ = m;4 + c i) 
and compare the relative merit of each choice through the 
calculated m% and d parameters. In Figure [7] we show how 
the measured ellipticites, after deconvolution, compare to 
the input ones for different choices of fixing /?. For clarity 
we have suppressed the error bars on all the curves barring 
those with the best and worst error bars. 

In all four cases we see that there is a strong dependence 
on the choice of f3. We see that in the eMERLIN simulations, 



the Massey & Refregier ( 2005 I optimisation algorithm per- 
forms badly in comparison to fixing n max and /3. In the 
LOFAR cases it performs a little better, but fixing the pa- 
rameters in question is still a marginal (~ 5% better on mi) 
improvement. The choice of /3 favoured in this approach 
is 0.2 x SExtractor FWHM; smaller j3 does not capture 
the information about shape at the edges of objects, while 
larger f3 smears out the detail in the object. We adopt 
/3 = 0.2 X FWHM for the rest of this work. 
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Figure 6. Whisker plots of point source ellipticities across the field for our four measurement sets, for both uncorrected (lower panels 
in each MS) and corrected ellipticities (upper panels in each MS). 
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Figure 7. Galaxy ellipticity, comparing input and measured ellipticities, for different choices of /3 for our four configurations. Black 
curves correspond to the shapelet optimisation algorithm of ( |Massey &: Rc frcgicr 2005]) , and the red, blue, green, yellow, brown and pink 
curves correspond to decompositions done with fixed n max — 10 and $ = 0.1, 0.2, 0.3, 0.4, 0.5 and 0.6 X SExtractor FWHM respectively. 
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3.5 uv Configurations and Impact on the PSF 

As stated in §2.5| we investigated the effect of regular and 
irregular sampling of the uv plane. For this investigation 
we have limited our analysis to the eMERLIN, 8 — +60° 
MS. Since the tracks in the uv plane are very similar for 
both samplings, the core of the PSF is unchanged; however 
the details of the uv sampling will alter the side lobes and 
structure of the PSF away from the main peak. 

In Figure [To] we show the PSFs generated by the simu- 
lation pipeline, resulting from the 2 different samplings dis- 
cussed above (top and middle panels), along with the dif- 
ference between the two (bottom panel). We see that the 
main lobe of both PSFs is very similar, and the difference 
at larger radii is relatively small. 

Having run simulations with the two different sam- 
plings, we matched the deconvolved shapes from both simu- 
lations to the original input catalogues, in order to directly 
compare the effect of the sampling difference on the resulting 
recovered shape distribution. In Figure [5] we show how the 
cllipticities of the galaxies measured using the regular and ir- 
regular sampling compare (left most column). We also show 
how the ellipticities in the two different samplings compare 
to the input ellipticities. 

We see that when we compare the input ellipticities to 
those measured using the regular and irregular samplings 
(middle and right most columns) the scatters are very sim- 
ilar, a\ = 0.049 and a\ = 0.049 for the irregular and 
a\ = 0.048 and erf = 0.049 for the regular sampling. This 
demonstrates that the different side lobe noise arising from 
the differing PSFs in the 2 sampling configurations does not 
significantly affect our ability to recover the input shape in- 
formation. 

In Figure [8] we also show how the matched ellipticities 
from the regular and irregular sampling compare against 
each other (left most column). We see that the two distri- 
butions are well correlated with correlation coefficients of 
pi = 0.886 ± 0.003 and p 2 = 0.891 ± 0.003; our shape mea- 
surement method recovers the shapes at a similar level in 
terms of the scatter in both cases. The scatter in the ellip- 
ticity difference between regular and irregular sampling, is 
around 40% smaller than the scatter in ellipticity difference 
between input and one of the output sets. 

3.6 Recovered Ellipticity Distributions 

We now turn to the main result of this work, which is to 
demonstrate how well we can recover the input ellipticities 
given the data reduction and imaging pipelines described 
above. We have described how we obtain measurements of 
input ellipticity e* and measured ellipticities et ; now we com- 
pare the two. 

In Figure [9] we show how our four MSs perform in this 
comparison (these are the blue curves in Figure [7|. For each 
of the four cases, we compute the Pearson's correlation co- 
efficient to quantify the correlation between the input and 
measured ellipticities, and also calculate the best fitting lin- 
ear model (ej — e| = m^e* + Cj) as well as a Xdof to assess 
the goodness of fit, all of which we summarise in Table |3] 

rrii 7^ is indicative of a calibration bias, which usually 
results from poor correction of factors that circularise im- 
ages, such as poor PSF correction. A non-zero value for a 
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suggests a systematic that induces some constant ellipticity 
so that even circular objects appear to have some ellipticity. 

Encouragingly we find that in all four of our test cases, 
we are able to recover tightly correlated input-to-output 
ellipticities, with the Pearson correlation coefficient being 
close to one in all cases (see Table [3). From Figure [9] we sec 
that for both LOFAR and eMERLIN we find a very close 
relationship between the original and measured ellipticities. 
The fact that there is essentially no Ci component to any 
of the simulations tells us that there is no constant induced 
ellipticity in our pipelines. There is evidence of a small cali- 
bration bias (i.e. non-zero m;) in all four sets of simulations, 
with the LOFAR ones faring slightly better than the eMER- 
LIN counterparts. 

We see that in both cases, the lower declination simu- 
lations provide a better recovered ellipticity (comparing mi 
values) over the DEC = +90° observations. The most accu- 
rately recovered ellipticity distribution is for LOFAR DEC 
= +60°, which we can see from Figure [3] has the fullest uv 
coverage, which should naturally lead to better quality im- 
ages. The worst performing simulation is eMERLIN DEC 
= +90°, which has the simplest uv coverage (c.f. Figure ^. 




Figure 10. PSFs obtained from the irregularly (top) and regu- 
larly (middle) sampled uv coverage. Also shown (bottom) is the 
difference between the two PSFs. 



3.7 Adding Noise 

The simulations above contain background fluctuations due 
to side-lobes, but do not contain measurement noise on the 
visibilities; we now consider the addition of this noise to 





Figure 8. Comparison of measured deconvolved and input ellipticities for an irregular and regularly sampled uv coverage. The leftmost 
column shows the ellipticity comparison between the 2 different samplings. The top (bottom) middle and rightmost panels shows the 
measured and input ellipticity comparison for the regular (irregular) sampling. Note this figure only shows the eMERLIN, 5 = +60° MS. 



Table 3. Table summarising ellipticity measurement fidelity, for all the reported telescope configurations. 



MS 


DEC (<5) 


NcLN 


Number Density 


n (arcmin 2 ) 


Pearson Correlation 


Best Fit Parameters 






(°) 


xlO 5 


Detected 


Useable 


Pi 


mi 


Ci 


eMERLIN 


90 


3 


128.0 


58.7 


0.720 ±0.008 
0.720 ±0.007 


-0.176 ±0.012 
-0.166 ± -2 x 10~ 4 


0.006 ± 0.004 
0.011 ± 0.004 


0.971 
1.400 


eMERLIN 


60 


3 


131.3 


70.1 


0.802 ± 0.005 
0.815 ±0.005 


-0.152 ±0.008 
-0.131 ± 0.008 


5 x 10" 5 ± 0.003 
0.003 ±0.003 


1.105 
0.852 


LOFAR 


90 


2 


121.4 


50.0 


0.687 ±0.009 
0.716 ±0.008 


-0.085 ±0.013 
-0.034 ± 0.012 


0.009 ± 0.005 
-0.001 ± 0.005 


0.843 
1.291 


LOFAR 


60 


2 


128.3 


59.3 


0.778 ±0.006 
0.816 ±0.005 


-0.068 ±0.008 
-0.075 ±0.007 


0.010 ±0.004 
0.003 ± 0.003 


0.961 
1.060 


LOFAR (N) 


60 


2 


127.1 


52.6 


0.759 ±0.007 
0.749 ±0.007 


-0.084 ±0.010 
-0.099 ±0.011 


0.008 ± 0.004 
0.002 ± 0.004 


1.270 
0.708 



the simulations. Within MeqTrees, the required input is the 
standard deviation of the Gaussian from which the noise 
for each visibility datum is drawn. For this test we restrict 
ourselves to the LOFAR DEC = ±60° measurement set. 



3.7.1 Assessing Noise Levels 

To assess the noise levels in the image plane corresponding 
to the noise input parameter in the MeqTrees software, we 
write the total noise nt in an image as 



2,2 
* ^sl ' 5 



(13) 
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Figure 9. Comparison of original and recovered ellipticities for our four configurations. 



where n s i is noise associated with side lobes (which is present 
even in the absence of measurement noise) and n m is the 
contribution arising from the corruption of the visibilities. 
We use the noise-free simulations described in the previous 
section to estimate the n s i term, and we run test simula- 
tions with varying noise in the visiblilities to estimate the 
n t terms. Hence we can estimate the n m terms using the 
equation above. 

We find that the relationship between the measured n m 
and the visibility noise is linear, and calibrate the chosen 
level of noise accordingly. We calculate the SNR of each 
object in the image by measuring the flux density in an 
aperture of 0.75" around its known position, and dividing 
through by the total noise n t integrated in the aperture. 

We choose the noise level such that our objects have 
SNR ~ 10; this is the level to which weak lensing measure- 
ments might plausibly be made. We found that in MeqTrees 
parameter units, a choice of 20 for the visibility noise cre- 
ates an image where the mean of the SNR distribution of 



the objects is (SNR) ~ 9.5, and this is the value that we 
adopted. 

With this choice of visibility noise we carry out the same 
procedure as before. We first determine how many CLEAN 
components we need to sufficiently deconvolve the beam, by 
running a similar test to the noise free case, the result of 
which is shown in Figure [TTj 

We find that there is a different behaviour to that seen 
in the noise free case (c.f. Figure [2]). We see that there is 
a reduction in the residuals as we increase the number of 
CLEAN components; at high numbers of components, both 
the point source and galaxy rms flattens out. In our noise 
free simulations we adopted a value of 1 x 10 5 and 2 x 10 
CLEAN components for point sources and galaxies respec- 
tively for this MS. We see that with the addition of noise, 
point sources require more CLEANing than before, while for 
the galaxies 2 x 10 iterations still seems sufficient. We do 
not find a minimum in the rms residuals for either the point 
sources or the galaxies, but since after Ncln = 2 x 10 there 
is only slow reduction, this is the value we have chosen. The 
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Figure 11. Mean rms residuals for varying numbers of Clark 
CLEAN components for the LOFAR DEC = +60° case. 



choice will be vindicated if there is good correlation between 
measured and input ellipticities. 

We now consider the noisy point source shape mea- 
surements. We perform a similar analysis for the shapelet 
modelling as in Section [3. 3| to find a good choice for the /3 
parameter. We decompose all the sources using the shapelet 
optimisation algorithm as well as using our fixed /3 and n max 
approach. Again, we use a fixed n max = 10 and use ft values 
of 1, 1.5 and 2 pixels. We show the corresponding results in 
Figure [12] 

Immediately we note that the addition of noise has 
broadened out the ellipticity distributions in contrast to 
those measured in Figure [5] As before, we find that fixing 
the j8 parameter seems to produce a cleaner measurement 
than the standard shapelet optimisation algorithm, and we 
see again that a choice of /? = 1.5 pixels produces the most 
tightly peaked histogram. We again adopt this value for f3 
and an n max — 10 in our modelling of the point sources and 
the PSF. 

We use the estimated PSF to deconvolve the noisy 
galaxies. From the previous section, we expect that the /3 pa- 
rameter choice for the galaxies is important, so we shapelet- 
decompose the galaxies with a variety of /3 choices. As before 
we find a strong /3 dependence when we examine how the 
measured ellipticities compare to the input ellipticities, as 
shown in Figure |13| again we only show the smallest and 
largest error bars for clarity. Similarly to the noise-free case, 
we find that fixing to 0.2 x the SExtractor FWHM es- 
timate yields the best results for m when we compare the 
measured and input ellipticities. We show in Figure [14] our 
final result of comparing the input and measured elliptic- 
ities in our noisy LOFAR DEC = +60° simulation, also 
summarised in Table [3] as LOFAR (N). 

We find that the addition of the visibility noise has 
degraded our slope measurement to ei — e* = (—0.084 ± 




0.10 



Figure 12. Normalised point source ellipticity histograms for 
different choices of /3, for the noisy LOFAR DEC = +60° case. 
Red curves correspond to the shapelet optimisation algorithm of 
( |Massey fc R cfrcgicr 2005), and the yellow, blue and green curves 
correspond to decompositions done with fixed n max = 10 and 
= 1, 1.5, 2 pixels respectively. The solid lines show e\ and the 
dashed show e 2 . 



0.010)ei + (0.008±0.004), and e 2 -£2 = (-0.099±0.011)e 2 + 
(0.002 ± 0.004). While this represents an increase in the cal- 
ibration bias over the noise-free case, the calibration bias 
remains at a modest 10% level. 



4 CONCLUSIONS AND FUTURE WORK 

In this paper, we have made an exploratory study of a ra- 
dio imaging pipeline suitable for studying weak gravitational 
lensing, and have tested a shear measurement pipeline which 
has been adapted for radio data. Using our simulations, we 
can obtain a better understanding of suitable shear mea- 
surement techniques to use with eMERLIN, LOFAR and 
ultimately the SKA. 

We have constructed a pipeline to simulate current and 
future weak lensing observations with radio interferometers. 



We were motivated by Patel et al. (20101 in which we at- 



tempted to measure a weak lensing signal using radio data 
from MERLIN and VLA. In that work we found systematic 
contamination in the data; indeed, a primary conclusion of 
that work was the need for a detailed study of the system- 
atics involved in trying to measure weak lensing using radio 
datasets. In the current work, we have been able to assess 
some of the possible systematics, and have demonstrated 
the reliability of our shape measurement method on realis- 
tic simulated radio images. 

Using the shapelets method we have created a set of 
images containing a collection of point sources and realistic 
galaxy shapes. The point sources were used to probe the 
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Figure 13. Ellipticity comparisons for different choices of /3 for the noisy LOFAR DEC = +60° MS. Black curves correspond to the 
shapelet optimisation algorithm of ( |Massey fe R cfrcgicr 2 005} , and the red, blue, green, yellow, brown, pink and purple curves correspond 
to decompositions done with fixed n max = 10 and /3 = 0.1,0.2,0.3,0.4,0.5,0.6 and 0.7 X SExtractor FWHM respectively. 
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Figure 14. Comparison of original and recovered cllipticities for the noisy LOFAR DEC = +60° MS. 



behaviour of the beam, while the galaxies were used to test 
how well one can recover known ellipticities in the presence 
of a radio data reduction and imaging pipeline. The 'true' 
images were run through our simulator to mimic observa- 
tions made with the eMERLIN and LOFAR arrays, at two 
different declinations, and with different numbers of CLEAN 
iterations. 



decomposition, including deconvolution of the PSF. As in 
our previous analysis, we have found that best results with 
shapelet decomposition were obtained by fixing the n max 
and /3 parameters. The galaxy simulations showed a very 
strong P dependance when we compared their true and 
measured ellipticities; we found that /3 = 0.2 x SExtractor 
FWHM gave the most faithful galaxy ellipticities. 



We measured ellipticities for galaxies via a shapelets 



Given our best shapelet models, we were then able to 
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compare the true and measured ellipticities in our cata- 
logues. We found highly correlated results, with all four MSs 
having Pearson correlation coefficients close to one. All MSs 
showed no evidence for an additive bias to the ellipticity 
measurements, and showed a modest (~ 10%) multiplica- 
tive bias. 

We added measurement noise to our simulations, fixing 
the visibility noise so that the resulting galaxies had a SNR 
distribution peaking at SNR ~ 10. In this case, we found 
similar results for the shapelet ellipticity measurements as 
with the simulations only containing side-lobe noise; the 
multiplicative bias for ellipticity measurement remains at 
the 10% level. 

These results are encouraging since they suggest that 
we can already recover well the true shape of sources from 
radio data, with well-motivated choices for the number of 
CLEANs, and the shapelet scale size. Given further studies 
of systematics effects, we hope to further reduce the calibra- 
tion bias. 

In this work we have demonstrated the feasibility of 
making weak lensing measurements in the presence of real- 
istic features of radio observations. Clearly, there are impor- 
tant extensions to what has been explored here: 

• In this study we have confined ourselves to small scale 
images, in order to CLEAN in a reasonably short time. IThe 
image size needs to be upscaled in order to probe position 
dependent PSF effects. 

• We have restricted ourselves to one CLEAN algorithm; 
the performance of a range of deconvolution techniques 
should be tested. 

• We have used a uniform weighting scheme; weight- 
ing this way should maximise the contribution from the 
longer baselines, resulting in better resolution images. Nat- 
ural weighting also exists, which give more sensitive images. 
Weak lensing is unique in that it requires both sensitiv- 
ity and high angular resolution images; there are weight- 
ing schemes that try to find a compromise between these 
two criteria (e.g. Briggs weighting, Briggs 1995). Assessing 
the shape measurement improvement/degradation between 
uniform, natural and intermediate weighting, and compar- 
ing this to any increase/decrease in source counts will be 
valuable. 

• In our construction of the MSs we have used particular 
frequency and time averaging configurations; these can dra- 
matically increase/decrease the size of the MS. Examining 
how different frequency and time averaging configurations 
affect shear estimates is a further important line of enquiry. 

• In this work we have used the shapelet method for 
shape measurement. The use of the shapelets method for 
radio data is well motivated, particularly by the Fourier in- 
variance of the shapelet basis functions and the possibility 
of shape measurement directly in the uv plane. However, 
in weak lensing, there are a considerable number of meth- 
ods for shape measurement, see for example Kitching ct al. 
( |2012[ > . Exploration of radio weak lensing should include a 
study of how these existing methods fare with radio image 
simulations. 

In this analysis we have demonstrated an approach to 
weak lensing measurements at radio wavelengths; encourag- 
ingly we have seen that our shape measurement methods 
and deconvolution techniques provide us with shape mea- 



surements that are only modestly biased. We have pointed 
out several open questions that still remain to be answered 
in this field, some of which this pipeline should be able to an- 
swer. These simulations can be extended to simulate actual 
weak lensing fields, where the input images contain realis- 
tic lensing shear, and we can assess how well we recover a 
cosmic signal. 

There are a range of new radio facilities that are/will be 
capable of producing data with which weak lensing measure- 
ments will be made. Understanding the relevant systemat- 
ics, and knowing how well shape-measurement methods per- 
form, is an important step towards using arrays such as a 
full Europe- wide LOFAR, and ultimately the SKA, for weak 
lensing. 
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